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ABSTRACT 



o 



Combining measurements taken using the Wilkinson Microwave Anisotropy Probe {WMAP) from 2001 to 2008 with measurements 
. taken using Planck from 2009 to 2010, we investigate the long-term flux density variability of extragalactic radio sources selected 

\^ \ from the Planck Early Release Compact Source Catalogue. The single-year, single-frequency WMAP maps are used to estimate 

yearly-averaged flux densities of the sources in the four WMAP bands: Ka (33 GHz), Q (41 GHz), V (61 GHz), and W (94 GHz). We 
identify 82, 67, 32, and 15 sources respectively as variable at greater than 99% confidence level in these four bands. The amplitudes of 
variation are comparable between bands, and are not correlated with either the flux densities or the spectral indices of the sources. The 
Q ' number counts of WMAP Ka-band sources are stable from year to year despite the fluctuation caused by individual source variability. 

5-H , Most of our sources show strong correlation in variability between bands. Almost all the sources that show variability are blazars. We 

■ have attempted to fit two simple, four-parameter models to the time-series of 32 sources showing correlated variability at multiple 

' frequencies - a long-term flaring model and a rotating-jet model. We find that 19 sources (60%) can be fit with the simple rotating-jet 

model, and ten of these also fit the simple long-term flaring model. The remaining 13 sources (40%) show more complex variability 
behaviour that is not consistent with either model. Extended radio galaxies in our sample show no sign of variability, as expected, 
with the exception of Pictor A for which we report evidence for a millimetre flare lasting between 2002 and 2010. 

> . 

. Key words, surveys: radio sources - radio continuum: galaxies - BL Lacertae objects: general - quasars: general - galaxies: individ- 

ual (Pic A) 

^ ' 

fS| '[ 1. Introduction detected by the Planck satellite during the first 1.6 sky surveys 

at all nine frequency channels. For some of the extragalac- 
Compact radio sources are the most significant contaminant j-^^jjo sources, the spectral shape at Planck frequencies 
^ ; at small angular scales in high-precision cosmic microwave suggests that they have been caught in a bright, flaring stage 
— .background (CMB) experiments such as the Wilkinson (Pianck Collaboration XIV 201 1). Knowledge of variability is 
Microwave Anisotropy Probe {WMAP, Bennett etal. 2003a) thus essential in order to correctly account for the contamination 
^ and the Planck' satellite (Planck Collaboration I 2011). In of the CMB measurements by radio sources. 
\^ , the frequency range 20-lOOGHz, which includes all the 
C3 . WMAP bands and the lower-frequency Planck channels. Most variable extragalactic sources are active galactic nu- 
" " the radio source population is dominated by compact, flat- clei (AGN), and the most extreme variability is seen in blazars. 
spectrum sources (Bennett et al. 2003b; Hinshaw et al. 2007; Relativistically beamed, non-thermal emission dominates the 
Wright et al. 2009; Planck Collaboration XIII 201 1), many of spectral energy distribution (SED) of blazars over other com- 
which are variable on time-scales ranging from days to years ponents such as dust emission from the galaxy or unbeamed 
(Planck Collaboration XIV 201 1 ; Planck Collaboration XV synchrotron emission from radio galaxy lobes. Relativistically 
2011). The Planck Eariy Release Compact Source Catalogue beamed sources are prone to variability for two reasons: (a) the 
(ERCSC, Planck Coflaboration VII 2011) lists bright sources emission is strongly Doppler boosted, so smafl changes of the 

Doppler factor caused by variations of the bulk Lorentz fac- 

1 „, ; ~ ~ ■ /m IN- ■ . r .u tor or the viewing angle can lead to strong variability; (b) in 

' Planck (http://www.esa.int/Planck) is a proiect of the ..... „ , . , . . V , . 

European Space Agency (ESA) with instruments provided by two sci- ^ relativistic jet, any small perturbation at the ongin of the jet 

entific consortia funded by ESA member states (in particular the lead can cause the formation of shocks (e.g., Begelman et al. 1984). 

countries France and Italy), with contributions from NASA (USA) and Variability caused by changes in Doppler boosting, also called 

telescope reflectors provided by a collaboration between ESA and a sci- geometric variability, changes the spectrum only slowly if at all 

entific consortium led and funded by Denmark. (c-g-, Begelman et al. 1980; Camenzind & Krockenberger 1992; 



X 



1 



X. Chen et al.: Long-term variability of extragalactic radio sources in the Planck ERCSC 



Steffen et al. 1995; Villata & Raiteii 1999). Flares produced by 
shocks in the plasma flow can lead to strongly frequency- 
dependent variability on shorter time scales, while also con- 
tribute to largely achromatic variability on longer time scales 
(e.g., Marscher & Gear 1985; Valtaoja et al. 1988). 

Studies of radio source variability in the radio-to-millimetre 
regime on time-scales up to several decades have been fa- 
cilitated by monitoring programs targeting medium or large 
source samples selected by known variability (Aller et al. 1996; 
Tornikoski et al. 1996; Aller etal. 1999; Stevens et al. 1994), 
radio spectrum (Valtaoja et al. 1992; Teraesranta et al. 1998; 
Nieppola et al. 2007), or gamma-ray brightness (Fuhrmann et al. 
2013). Discussions of blazar variability have mostly been fo- 
cused on understanding the characteristics of flares, which 
are observed on time scales ranging from one day or less 
(Wagner et al. 1995) to weeks or months (e.g., Angelakis et al. 
2012) to several years (Hovatta et al. 2009), and their inter- 
pretation in the framework of jet models (e.g., Valtaoja et al. 
1999; Liihteenmaki et al. 1999; Lahteenmaki & Valtaoja 1999; 
Tiirler et al. 2000; Ciaramella et al. 2004; Hovatta et al. 2009; 
Angelakis et al. 2012). While it is generally accepted that 
shock-produced flares dominate the variability on very short 
time scales (days to weeks), geometric variability should have 
an effect on longer time scales because rotation or preces- 
sion, as suggested by VLBI observations of many jets (e.g., 
Kellermann et al. 2004), would inevitably lead to changes in 
Doppler boosting of the radiative zones in the jet. Indeed, 
Bach et al. (2006) have modelled optical and radio time-series 
of blazars observed in the monitoring campaigns of the Whole 
Earth Blazar Telescope (WEBT, Villata et al. 2002, 2004) in the 
framework of inhomogeneous, rotating helical jets (Villata et al. 
1998; Villata & Raiteri 1999; Ostorero et al. 2004), and have 
suggested that all variability down to time scales of ~ 100 days 
can be explained by geometric effects alone. 

In this paper, we report results on long-term flux density vari- 
ations of a sample of sources selected from the Planck ERCSC. 
We primarily use data collected by WMAP and Planck. WMAP 
started its observations at the second Sun-Earth Lagrange point 
(L2) in August 2001, and ended the collection of science data 
on 19 August 2010. Planck started its observation at L2 on 13 
August 2009 and is still taking data. Both telescopes take about 
6 months to complete a single sky survey. Since the lower four 
frequency bands of Planck overlap with the upper four frequency 
bands of WMAP, and both telescopes have the powerful ability 
to make near-simultaneous measurements at multiple frequency 
bands, combining the observations from the two experiments 
across an interval of nearly ten years provides the best oppor- 
tunity to track the long-term variations of radio sources at fre- 
quencies between 30 and 100 GHz. Table 1 lists the central fre- 
quencies and beam-widths of the WMAP and Planck bands that 
are of interest here. 

We describe the selection of our sample from the Planck 
ERCSC in Sect. 2. In Sect. 3 we introduce the WMAP data used 
in our analysis and the method of estimating source flux densities 
from the WMAP maps. Source variability is quantified and stud- 
ied statistically in Sect. 4. In Sect. 5, we analyse the long-term 
averaged variability patterns of strongly variable blazars in our 
sample. We then discuss the interesting case of long-term vari- 
ability in unbeamed radio sources in Sect. 6. The results are put 
into perspective with other studies in Sect. 7, and we conclude 
in Sect. 8. 



2. Sample Selection and Planck Data 

The Planck Early Release Compact Source Catalogue 
(Planck Collaboration VII 2011) is a catalogue of highly 
reliable sources (with > 90% cumulative reliability), both 
Galactic and extragalactic, detected over the entire sky in the 
first 1.6 Planck sky surveys. This catalogue provides separate 
lists of sources detected in each of the nine Planck frequency 
channels (30 to 857 GHz). We restricted our sample selection to 
the Planck 30, 44, 70, and 100 GHz channels that overlap with 
the WMAP Ka to W bands, and excluded sources close to the 
Galactic plane (|^| < 5°). We first collated the lists at these four 
frequencies and constructed a "band-merged" catalogue listing 
the sources detected in all four channels. The band-merging 
process was based on positions only. No brightness information 
was used, to avoid biasing towards or against any class of 
sources. The source list at each frequency was used as the 
seed catalogue and counterparts were sought in the candidate 
catalogues, i.e., the source lists at the other frequencies. The 
matching radius was always defined as half of the full width at 
half maximum (FWHM) of whichever band (seed or candidate) 
has lower resolution. A confusion check was then performed 
to ensure a one-to-one reciprocal relation between matches 
across all four bands. Due to the change in resolution, there 
was one case in which the 44 GHz source could be matched 
with two sources at 100 GHz; we excluded this source. We then 
cross-checked this list with both the NASA/IPAC Extragalactic 
Database (NED-) and the SIMBAD Astronomical Database^^ to 
remove any Galactic objects (Hii regions, planetary nebulae, 
supernova remnants, etc.). This process yielded a final list of 
198 extragalactic sources, including 135 flat-spectrum radio 
quasars (FSRQs), 26 BL Lacertae objects (BLLs), 28 non-blazar 
active galaxies (AGNs, including many extended radio galax- 
ies), 1 starburst galaxy (SBG), 1 normal galaxy (GAL), and 7 
radio sources of uncertain type. The classification of the active 
galaxies was taken from the Candidate Gamma-Ray Blazar 
Survey (CGRaBS, Healey et al. 2008), the third edition of the 
Roma-BZCAT blazar catalogue (Massaro et al. 2011), and the 
twelfth edition of the Catalogue of Quasars and Active Nuclei 
(Veron-Cetty & Veron 2006). Fig. 1 shows the distribution of 
the sources on the sky. We emphasize that there is no guarantee 
that the sample is complete. In addition, requiring that a source 
be seen in all the Planck bands between 30 and 100 GHz will 
inevitably favour flat-spectrum radio sources. 

The calibration of the Planck 30-100 GHz maps is 
currently based on the CMB dipole (Zacchei et al. 2011; 
Planck HFI Core Team 2011). Since the ERCSC flux densities 
were calculated directly from the maps and are therefore only 
correct if the source has a CMB spectrum, colour corrections are 
needed for sources with different spectra. To determine the spec- 
tra, we looked for the counterparts of our sources in the higher- 
frequency (143-857 GHz) ERCSC source Usts. We then fit the 
SED of each source over all the frequencies at which it was de- 
tected, using three simple models, when there were enough data 
points for fitting: 

i) single power law 

log 5 = P0+ pilogv; (1) 

ii) quadratic fit 

log 5 = /7o + f 1 log V -H p2(\og vf ; (2) 

^ http://ned.ipac.caltech.edu/ 

^ http : // simbad . u- strasbg . fr/ simbad/ . 
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Table 1. Characteristics of relevant WMAP and Planck bands in this study. 





WMAP 






Planck 




Band 


Center Frequency [GHz] 


Beam FWHM [arcmin] 


Band 


Center Frequency [GHz] 


Beam FWHM [arcmin] 


K 


22.8 


49.2 








Ka 


33.0 


37.2 


30 


28.5 


32.65 


Q 


40.7 


29.4 


44 


44.1 


27.92 


V 


60.8 


19.8 


70 


70.3 


13.01 


w 


93.5 


12.6 


100 


100.0 


9.37 



Notes. Values for WMAP are from Bennett et al. (2003a) and those for Planck are from Planck Collaboration I (201 1). 



+ 90 




-90 



Fig. 1. Map in Galactic coordinates showing the locations of the 198 extragalactic sources selected from the Planck ERCSC. The 
north ecliptic pole (NEP) and the south ecliptic pole (SEP) are indicated by crosses. Source densities are higher near the ecliptic 
poles, where Planck has longer integration time. 



iii) double power law 

S = Po(v/vor + Piiy/voY' ■ (3) 

Here /7;(i=0, 1,2,3) denote parameters that are different in 
each model. Since the Planck 100 GHz data can be signif- 
icantly contaminated by the CO (7=1— *0) line at 115 GHz 
(Planck HFI Core Team 201 1), we did not include the 100 GHz 
flux densities in the SED fitting. A model fit was accepted when 
the;^'^ values indicated that the model was compatible with 95% 
of the data. In cases when more than one model was accept- 
able, we always selected the model with fewest parameters un- 
less an additional parameter significantly improved the value of 
the reduced x~ (as described in Bevington & Robinson 2003). 
A few examples of the source SEDs and our fitted models are 
given in Fig. 2. For sources with a good fit (144 out of the 
198 sources), we used the model to estimate the spectral index 
a (S cc v") at the Planck central frequencies. For sources that 
could not be reasonably fit by any of the models, flux densities 
in the adjacent bands were used to estimate the spectral indices. 



We then applied the corresponding colour corrections given in 
Zacchei et al. (201 1) and Planck HFI Core Team (201 1). These 
range from almost no change to 5.7%, depending on frequency 
and source spectral shape. The absolute calibration errors are at 
the 1% level for the Planck 30-70 GHz channels, and 2% level 
for the 100 GHz channel. We conservatively assumed an overall 
calibration uncertainty of 3% and added it in quadrature with the 
ERCSC flux density errors prior to modelling the spectra. 



3. WMAP Data 

The WMAP individual-year maps from 2001 to 2008 and the 
flux measurements of our sources from these maps constitute the 
main dataset of this paper, which is described in this section. To 
support our analysis, we have also used some ground-based data 
from the Effelsberg 100-m telescope, the IRAM 30-m telescope, 
and the Metsahovi 14-m telescope. We defer the description of 
these data to the sections where they are used. 
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Fig. 2. Examples of source SEDs and fitted models (grey solid 
line). The reduced values are indicated. 



3.1. WMAP maps and photometry 

The high-resolution (A^side = 1024 in the HEALPix scheme, 
Gorski et al. 2005) WMAP single-year, single diff'erencing- 
assembly (DA) maps^ from 2001 to 2008 are the basis of this 
study. Since the ERCSC was derived from single-frequency 
Planck maps, we performed a weighted, pixel-by-pixel, mean 
of the single DA maps from each year to generate single-year 
maps at each WMAP frequency. The number of observations per 
pixel was used as the weight. The average temperature outside 
the WMAP seven-year KQ85 mask (Gold et al. 201 1) was sub- 
tracted from the single DA maps before the weighted mean was 
calculated. 

For consistency with the ERCSC flux density measurements, 
we used aperture photometry to obtain flux densities for the 
PZfl«cA:-detected sources from the WMAP maps. The source posi- 
tion was taken from the ERCSC 100 GHz catalogue as this band 
has the highest resolution of the four bands considered here. 
We measured the flux density in a circular aperture with radius 
equal to 1.1 times the nominal FWHM of the WMAP beam at 
each frequency, and subtracted an estimate of the local back- 
ground, computed as the median intensity within a concentric 
annulus. The size of the aperture is chosen to include > 95% 
of the source power in the K- to V-bands, which was calculated 



All the maps are available from NASA's Legacy Archive 
for Microwave Background Data Analysis (LAMBDA) at 
http : //lambda . gsf c. nasa . gov/ . 



from the WMAP seven-year beam radial profile^. Due to a signif- 
icant shoulder in the W-band beam profile at 0°2-0?5 (Page et al. 
2003), this aperture only includes ~ 85% of the source flux in 
the W-band. To have > 95% of the source power in the aperture, 
we would have to adopt an aperture radius of 2xFWHM in this 
band, which would then include too much noise from the back- 
ground. We therefore decided to stay with an aperture radius of 
1 . 1 xFWHM for the W-band, and correct for the flux loss using 
Monte Carlo simulations (see the next paragraph). Similarly, we 
calculated at each band the radial distance from the center of the 
beam that includes 97.5% of the source power (~1.3, 1.3, 1.4, 
1.7, or 2.4 times the corresponding beam FWHM at K, Ka, Q, 
V, or W band), and used this for the inner radius of the annulus 
to avoid over-subtraction of the source. We then chose the outer 
radius to make the area of the annulus approximately the same 
as the aperture. 

For studying variability, it is important to have unbiased flux 
density measurements and reliable estimates of their uncertain- 
ties. Since we do not know whether the aperture is perfectly cen- 
tered on each source and we do not have precise estimates of the 
background, we estimated aperture corrections using a Monte 
Carlo method. We injected fake point sources into the WMAP 
maps at random positions and compared the extracted flux densi- 
ties with the input. The flux densities of the artificial sources had 
a flat dlogN/dlogS distribution (equal numbers of sources per 
decade of 5) in the range 100 mJy < S < 200 Jy. The sources 
were first injected into single pixels in an empty A^side = 4096 
map, which was then smoothed with the WMAP beam profile ap- 
propriate for the chosen frequency and degraded to A^side = 1024 
before being added to the WMAP single-year maps at that fre- 
quency. We then performed the same aperture photometry at the 
input source locations, excluding those at \b\ < 5° and those 
within 2° of other fake sources or sources in the WMAP seven- 
year source catalogue (Gold et al. 201 1). To minimize collisions 
of fake sources, we injected 1000 sources in each map and it- 
erated five times to get a sample of ~ 5000 fake measurements. 
We used the distribution of the ratio of extracted to input flux 
density to estimate a bias correction factor, which was typically 
~ 10-20% except at the W-band. 

We estimated the uncertainties in the flux densities obtained 
from the aperture photometry by error propagation in the for- 
mula for the aperture flux density, assuming white uncorrelated 
noise. This assumption is justified for instrumental noise, but not 
for noise in the background. We therefore used Monte Carlo sim- 
ulations to test our error estimates, by injecting 1000 sources 
with the same flux density (from 1 Jy to 20 Jy in 1 Jy steps) in 
each map, and using the dispersion of the recovered flux densi- 
ties as an indicator of the "true" uncertainty. We found that the 
flux density uncertainties from the white noise model were al- 
ways less than those returned from Monte Carlo simulations, as 
expected. However, when repeating the same process on CMB- 
subtracted maps (subtracting the WMAP seven-year internal lin- 
ear combination CMB map at \b\ > 30°), the simulation results 
are in good agreement with the white-noise estimates, showing 
that the correlated noise is mostly CMB. Since CMB noise is 
constant from year to year and therefore will not affect the vari- 
ability analysis, we decided to use the flux density uncertainty 
returned directly from the white-noise model in our study. 

The absolute calibration of WMAP data is based on the CMB 
monopole temperature and the velocity-dependent dipole result- 
ing from WMAP's orbit about the solar system barycentre, and 
has an absolute calibration uncertainty of 0.2% (Jarosik et al. 



' Also available from NASAs LAMBDA site. 
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201 1). The aperture correction error measured from simulations 
varied between 0.4% and 1.3% in the K to W bands. We there- 
fore assumed an overall 2% error from aperture correction and 
calibration, and added it in quadrature with the flux density un- 
certainty to constitute the final error term for the flux density of 
each source. 

We fit the SED of each source at each year with a single 
power law (Eq. 1), and applied colour corrections to the central 
frequencies, following the recipe of Jarosik et al. (2003). The 
WMAP single-year maps are generally much noisier than the 
Planck maps, especially at the V and W bands where even the 
seven-year co-added maps are less sensitive than the one-year 
Planck maps (see Planck Collaboration VII 201 1, Fig. 5). As a 
result, the flux density estimates can sometimes be negative for 
faint sources at certain frequencies in some years. In the SED 
fitting, we only used the good quality (defined as 5 > 3cr) flux 
density measurements. A 3cr threshold (~ 99.7% probability that 
the signal is truly positive assuming the background is Gaussian) 
is adopted and considered sufficient in our case as we look for a 
positive signal at the location of a known source. We do, how- 
ever, caution that in non-Gaussian backgrounds a 3cr threshold 
raises the risk of deriving incorrect flux density and error es- 
timates. Although the WMAP K-band (23 GHz) flux densities 
were not used in the variability study because there is no corre- 
sponding channel in Planck, we included them in the SED fitting 
to obtain more accurate spectral indices. 

The WMAP single-year flux densities for our sources, to- 
gether with the Planck ERCSC flux densities extrapolated to the 
WMAP centre frequencies (using the SED models described in 
Sect. 2), are listed in Table 2. These flux density values are used 
to construct light curves for all the sources in our sample. The 
full set of light curves is available in the online supplementary 
material, arranged in ascending order of Right Ascension. Only 
light curves that have four or more good quality measurements 
in the seven years of WMAP data are plotted (see Sect. 4 for 
detailed discussion). For sources that are specifically discussed 
in the paper, individual light curves are presented where they 
are discussed. We show all the flux density measurements with 
S > 2cr, with dash-dot lines indicating 3cr noise level at each 
frequency. Since the flux measurements here are generally aver- 
ages over a fraction of a year (see Sect. 3.2), all the data points 
are plotted at a nominal epoch in February of each year, as both 
WMAP and Planck single year observations start in mid- August. 
The effective epoch of each observation may differ by a few 
months from the nominal epoch. 

3.2. Difference between WMAP and Planck flux density 
measurements 

WMAP and Planck employ different scanning strategies, and 
when comparing flux density measurements from the two tele- 
scopes it is important to understand how they are affected by the 
scanning strategy. 

The WMAP satellite is a differential experiment, and its opti- 
cal system consists of two back-to-back telescopes that measure 
the difference in temperature between two points separated by 
141° on the sky. The satellite spins around its symmetry axis 
once every 2.2 min while the spin axis precesses at 22.5° h"' 
about the Sxxn-WMAP line. Since each telescope line of sight 
is ~ 70° off the symmetry axis, the combined motion of spin and 
precession causes the observing beams to fill an annulus cen- 
tered on the local solar vector, with inner and outer radii of ~ 48° 



and ~ 93°, respectively. Although this complex scanning strat- 
egy provides many repeated observations ( WMAP observes more 
than 30% of the sky each day), the coverage is not uniform. For 
example, in one year, a source on the ecliptic plane can be seen 
continuously for 1.5 months, then off for 3 months, then on for 
another 1.5 months, then off for 6 months; while a source at 48° 
ecliptic latitude or higher can be seen for 6 months straight, and 
then not seen for another 6 months. Sources close to the ecliptic 
poles are seen continuously throughout the year Therefore, flux 
densities estimated from a single-year map are averages over a 
different fraction of the year for each source. 

The Planck satellite spins around its axis at a rate of one 
rotation per minute. The focal plane is oriented at an angle of 
85° to the spin axis, which tracks the direction of the Sun at 
~ 1° per day (Planck Collaboration I 2011). This scan pattern 
also results in inhomogeneous integration time, and areas near 
the ecliptic poles are observed with greater depth than the rest of 
the sky. Most sources are seen for a few days every six months. 
Since the data included in the ERCSC amount to 1.6 full-sky 
surveys (Planck Collaboration VII 2011), some of the sources 
in our sample have been covered twice with a time separation of 
~ 6 months. Therefore, with the exception of sources close to the 
ecliptic poles, the flux densities given in the ERCSC are mostly 
from a single snapshot or an average of two passes, rather than 
an average over a substantial fraction of a year 

Because of these differences, we treat the WMAP single-year 
flux densities and the ERCSC flux densities differently in our 
analysis. To avoid systematic errors, we include only the WMAP 
flux densities in the statistical studies of the source variability. 
The ERCSC flux densities, extrapolated to the WMAP frequen- 
cies, while helpful in confirming the long-term trend for sources 
exhibiting long-term variability, are more sensitive to short-term 
variability (e.g., flares). An example is given in Fig. 3, where the 
top panel shows that the highly variable source 3C273 main- 
tains a steep spectrum in the WMAP year 1 to year 7 observa- 
tions, but its spectrum suddenly becomes flat in the ERCSC. 
This is because this source was undergoing a strong flare dur- 
ing the time Planck observed it (1-10 January 2010). Ground- 
based observations with the Effelsberg 100-m and IRAM 30-m 
telescopes (obtained within the framework of the F-GAMMA 
program; Fuhrmann et al. 2013) show that this source has un- 
dergone a few flares in the years 2007-2008 while generally 
preserving a steep spectrum (Fig. 3, bottom panel). The January 
2010 observation at the Effelsberg telescope was obtained near- 
simultaneously with the Planck data, and confirms the flaring 
activity captured by Planck. 

4. Flux Density Variability 

In this section we quantify and discuss the variability properties 
of the sources in our sample. We first define and discuss dif- 
ferent measures of variability in accordance with those used in 
the literature, then investigate possible effects of variability on 
source number counts, and finally discuss correlation between 
variability found in different bands. In all steps of our analysis, 
we define incremental quality cuts to select sources which are 
included in the analysis. These cuts are chosen in order to max- 
imize the number of sources included while keeping the results 
robust. We discuss qualitative statistical trends found in our sam- 
ple, but emphasize that no quantitative significance levels can be 
assigned to them owing to the incompleteness of the sample. 

* Seethe VKMAP mission site (http: //wmap . gsfc .nasa. gov/) for 
more details of scanning strategy. 
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Table 2. Flux densities at WMAP Ka, Q, V, and W band for the 198 extragalactic radio sources selected from the ERCSC (Example, full table can be retrieved online). 



Name" 


GLON 


GLAT 


Type 


Band 


WMAPl 


WMAP2 


WMAP3 


WMAP4 


WMAP5 


WMAP6 


WMAP7 


ERCSC* 




Notes'* 




[deg] 


[deg] 






[Jy] 


[Jy] 


[Jy] 


[Jy] 


[Jy] 


[Jy] 


[Jy] 


[Jy] 






J0006-0623 


93.503 


-66.636 


BLL 


Ka 


1.44+0.33 


1 88+0 30 

J. • u _!_ yj u ^yj 


1.86+0.29 


1.54+0.29 


1.26+0.29 


1 05+0 30 

J. .yjy _!_ vy . y yj 


0.77+0.30 


1.46+0.08 


n.24 












o 


1.69+0.37 


1.32+0.34 


2.31+0.34 


1.74+0.35 


1 50+0 35 

± 1 .yyj-i—yjt .y y 


0.67+0.41 


1.20+0.36 


1.43+0.07 


10.96 












V 


1.72+0.58 


1.23+0.54 


1.19+0.63 


1.14+0.56 


1.49+0.59 


0.90+0.56 


1.27+0.63 


1.35+0.09 














w 


-0.04+0.93 


0.45+0.84 


0.48+0.91 


1.17+0.85 


3.95+1.09 


-0.37+0.72 


1.25+0.84 


1.22+0.08 






JOOlO+1058 


106.981 


-50.622 


AGN 


Ka 


0.65+0.29 


1.31+0.30 


1.60+0.32 


2 09+0 30 


90+0 29 

yj , yyj _!_ yj . jl— y 


92+0 30 

yj . y A^y—yj , y yj 


0.20+0.30 


2.51+0.10 


29.68 


V 


(Mrk 1501) 








o 


0.38+0.37 


0.70+0.37 


1.54+0.34 


2.11+0.37 


1.04+0.36 


0.53+0.37 


1.75+0.36 


2.68+0.16 














V 


-0.13+0.63 


0.75+0.56 


1.73+0.63 


0.75+0.60 


0.90+0.61 


0.54+0.64 


1.93+0.60 


2.27+0.17 




( 










w 


69+0 99 


1.07+0.97 


2.20+1.12 


3 54+0 75 


2.01+0.93 


-1.1 1+0.91 


0.98+1.01 


1.89+0.10 






J0050-5738 


303 285 


-59.481 


FSRQ 


Ka 


97+0 23 


0.94+0.22 


92+0 23 


1.61+0.22 


1 38+0 23 


1.64+0.22 


1 55+0 23 

I. ,y y _!_ yj » z. y 


1 33+0 09 

L ^y y _!_ yj ^yj y 


18.85 


V i 


(0047-579) 








o 


1.30+0.27 


1.21+0.25 


1.77+0.26 


1.62+0.27 


1.27+0.26 


1.14+0.27 


1.45+0.27 


1.36+0.18 


-9.43 












V 


1.46+0.42 


1.03+0.45 


1.21+0.45 


1.40+0.46 


0.57+0.43 


0.42+0.44 


2.08+0.45 


1.25+0.15 














w 


0.24+0.72 


1.53+0.73 


0.94+0.57 


1.71+0.81 


1.24+0.69 


52+0 68 

yj 1 y A^y—yj , yjyj 


0.04+0.68 


81+0 08 

vy - vj ± _!_ vy - vy vj 




d 


J005 1-0650 


122.730 


-69.711 


FSRQ 


Ka 


1 53+0 30 


1.11+0.30 


1.44+0.31 


1 36+0 30 

J. • ^1 yj _!_ yj u ^ yj 


1.16+0.29 


1.24+0.30 


2.01+0.30 


1.38+0.08 


5.43 




("0048-071") 








o 


1 70+0 38 


0.77+0.36 


1.06+0.37 


1.17+0.38 


-0.41+0.38 


76+0 38 

yj , 1 yj _!_ yj . y yj 


1 38+0 35 

1 . %y u _!_ vy . y y 


1.49+0.09 














V 


2 25+0 60 


1 25+0 63 


1.58+0.57 


0.44+0.58 


2.32+0.59 


2.62+0.60 


3 78+0 68 

y. 1 yj _!_ yj .yjyj 


1.55+0.11 


12.76 












W 


2.09+0.90 


0.66+0.94 


1.94+0.90 


2.70+0.83 


1.94+0.85 


-0.89+0.92 


2.05+1.16 


1.39+0.09 




1 


JO 106-4034 


290.673 


-76.187 


FSRQ 


Ka 


2.37+0.22 


1.42+0.22 


1.90+0.22 


2.64+0.22 


3.77+0.22 


3.51+0.22 


3.40+0.22 


1 30+0 07 

X • y yj _!_ yj • \y / 


31.75 


"V i 


("0104-408") 








o 


2.34+0.28 


1.67+0.28 


2.10+0.29 


2.60+0.28 


3 55+0 28 


3 39+0 28 

y ,y y _ !_ vy . 


3.17+0.29 


1 33+0 08 

± ^y y _!_ \y • \y vj 


24.00 


"V 










V 


1.47+0.45 


1.78+0.47 


2.24+0.49 


1.26+0.49 


3.39+0.44 


3 50+0 46 

y . yyjy—yj.^yj 


2.47+0.45 


1.23+0.11 


30.58 


"V < 










w 


3 27+0 60 

/ _!_ \J u \J\J 


1.57+0.62 


1.57+0.75 


3.22+0.67 


2.74+0.63 


2.67+0.73 


2.51+0.76 


97+0 08 

yj - y 1 _!_vy-vyvj 


-15.32 




J0108+0135 


131.817 


-60.985 


FSRQ 


Ka 


2 33+0 30 


3 32+0 31 


2.43+0.30 


1.12+0.31 


1.45+0.32 


1.05+0.30 


1.56+0.31 


3.13+0.10 


40.77 


V i 


(0106+013) 








o 


2.05+0.36 


3 06+0 37 


1.39+0.36 


1.87+0.36 


0.21+0.34 


0.47+0.37 


09+0 39 

yj ,yj y _!_ vy » ^ y 


2.85+0.08 


37.99 


V i 










V 


2.13+0.56 


2 36+0 62 

. ^yj 


2.31+0.61 


1.64+0.63 


46+0 63 

yj . ^vj_!_ vy . yj ^ 


0.67+0.62 


1 68+0 57 


2 38+0 06 

^^yyj _!_ vy - vy vy 














W 


-1.17+1.15 


1 89+1 06 


-1.34+1.17 


-1.63+1.20 


-0.34+0.73 


1.49+0.91 


0.07+1.20 


1.97+0.05 




1 


J0121+1 149 


134.577 


-50.354 


FSRQ 


Ka 


0.94+0.32 


0.82+0.31 


66+0 32 


0.99+0.31 


39+0 33 

yj ,y y _!_ vy . y y 


1.77+0.31 


1.81+0.31 


1.75+0.08 




1 

1 


^0119+115") 








o 


38+0 39 


0.92+0.36 


1.03+0.36 


0.83+0.37 


81+0 35 

yj ,yj \. _!_ vy . y y 


1.57+0.35 


2.65+0.37 


1.56+0.07 




1 
1 










V 


74+0 63 


-0 16+0 60 


95+0 66 

\j ~ y .y _!_ yj .yjyj 


63+0 63 

yj M\j.y _!_ w - vj^ 


1.14+0.59 


1 03+0 59 

\. .yj y _!_ vy . y y 


30+0 58 

vy . ^ vy — 1— vy . y u 


1.26+0.05 














w 


65+0 81 


-0.79+0.94 


15+0 89 

V/ - J. ^ _!_ yj y 


-1 63+0 89 

A. .yj^ _!_ vj. yj y 


0.27+0.77 


-0 40+0 93 

yj . 1 \_/ — !— V-/ » y y 


0.84+1.07 


1 00+0 05 

X - vy vy _!_ vy .yj y 






10132-1655 


168.145 


-76.020 


FSRQ 


Ka 


2.10+0.28 


2 56+0 30 


2 38+0 30 


2.47+0.30 


2.18+0.29 


1.82+0.29 


1.87+0.29 


1.71+0.08 


-3.11 




(■ono-171") 








o 


1 66+0 33 

J. . vJ V_/ _]_ yj t^,y 


2 09+0 33 


2 20+0 36 

^ • w _]_ yj t ^yj 


1 76+0 36 

L t 1 yj —1- yj t ^yj 


2 16+0 35 


1.54+0.34 


1.29+0.34 


1 53+0 06 

X • .y ^ _!_ vy . vy vy 


0.73 












V 


95+0 51 


2.71+0.55 


2.28+0.54 


2.77+0.56 


1 60+0 58 

X 1 yjyj -L- yj . y yj 


2.21+0.58 


1.24+0.52 


1.24+0.04 


-7.05 












w 


1 56+0 84 

± . ^ yj _!_ yj u u ~ 


4.50+0.77 


38+0 78 


1.54+0.98 


0.73+1.01 


1.49+0.75 


2 62+0 89 


99+0 04 

yj ~ y y _!_ vy . yii 






10136+4751 


130 788 


-14.316 


FSRQ 


Ka 


4 s^+o 28 


4 86+0 28 


4 61+0 28 

^ . yj I _1_ yJ t^yj 


4.13+0.29 


3 50+0 28 

y . yyj j-yj.-^yt 


3 66+0 29 


4 32+0 29 


4 68+0 12 

^ • VJO -I- Vy • X ^ 


9.59 


V ' 


(0133+476) 








Q 


5.02+0.37 


5.42+0.38 


4.86+0.36 


4.24+0.39 


2.42+0.35 


4.10+0.36 


4.01+0.36 


4.42+0.11 


21.99 


V 










V 


4.15+0.58 


6.32+0.56 


3.89+0.57 


3.39+0.60 


2.98+0.56 


3.44+0.61 


3.33+0.63 


3.92+0.12 


25.29 


V i 










w 


3.15+0.90 


2.49+0.69 


1.73+0.85 


1.80+0.93 


0.07+0.95 


0.66+1.09 


3.56+0.90 


3.40+0.11 




1 
1 


J0137-2431 


201.403 


-79.286 


FSRQ 


Ka 


0.53+0.28 


1.44+0.26 


1.43+0.28 


1.07+0.24 


1.23+0.27 


1.05+0.27 


0.57+0.27 


2.45+0.09 


-8.72 




(0135-247) 








Q 


0.82+0.33 


1.72+0.32 


1.37+0.31 


1.64+0.31 


1.27+0.33 


1.20+0.33 


1.41+0.33 


2.49+0.15 


-14.46 












V 


0.80+0.53 


1.12+0.50 


3.21+0.56 


1.72+0.52 


1.52+0.50 


0.98+0.52 


1.69+0.54 


2.73+0.14 


14.00 












w 


-1.96+0.81 


2.39+0.82 


2.73+1.01 


1.47+0.81 


1.30+0.81 


-1.11+0.75 


-0.23+0.83 


2.25+0.09 







Notes. <'''TheJ2000 name for each source is given. For easy reference, we have also included in parentheses the commonly known name or its B 1950 name. * ' These are the ERCSC flux densities 
extrapolated to the center frequencies of WMAP Ka, Q, V, and W-band. This field is left blank when a timeline is not included in the variability analysis. Negative values indicate that the flux 
density fluctuation level is less than the noise level (see Sect. 4.1). "V" denotes strongly variable sources (>99% confidence), and "v" denotes moderately variable sources (>95% confidence) at 
a given frequency. 
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Fig. 3. Top: WMAP/Planck light curves of 3C273. The WMAP 
flux densities are shown by filled symbols and the Planck flux 
densities are shown by open symbols. Bottom: spectra of 3C 273 
from Planck, Eff'elsberg, and IRAM observations at several 
epochs. The 201001 Eff'elsberg 100-m observation is nearly si- 
multaneous with the Planck observation, and confirms the flare 
seen by Planck. 

4.1. Quantifying Variability 

We define the fractional variability index, Vims, by 



where 



N- 



1 ^ <5>„ 



'-Y 



(4) 



(5) 



is the mean of the individual flux measurements 5,, weighted 
by their individual inverse variances corresponding to the mea- 
surement error cr, , and 



2,1/^ 



(6) 



is the average error of the measurements contributing to the 
weighted mean. The quantity is defined as 



1=1 



0"7 



(7) 



where 5 = {S}^, andx^/{N- 1) is the reduced of the sample 
for N - I degrees of freedom (assuming Gaussian errors). The 
term X^js gives the variability amplitude in units of the average 
measurement error As explained in Appendix A.l, our Vrms is 
equivalent to the quantity used by Sadler et al. (2006), but gen- 
eralized for the case of small sample size and unequal measure- 
ment errors. 

As defined in Eq. 4, Vrms naturally connects the variability 
amplitude of a source with the statistical significance of its be- 
ing variable (see Appendix A. 2). Following Bolton et al. (2006) 
and Franzen et al. (2009), is used to test the null hypothesis 
that each source is non-variable (i.e., having constant flux den- 
sity) at each frequency. We classify a source as strongly vari- 
able when the value indicates that the probability of the null 
hypothesis is less than 1% (i.e., source is variable at >99% 
confidence level). For N - 1 (corresponding to seven years of 
WMAP measurements), we have - 1 = 6 degrees of free- 
dom, which requires > 16.8 or Vrnis > 1.34 (cr)„ / (5 )w in 
order to claim a source strongly variable. Additionally, we use 
Vmis > <cr)^/(5)„ (or XiTOs > 1), corresponding to a confi- 
dence level of ~95%, to classify a source as moderately vari- 
able. Other sources are not considered to be significantly vari- 
able. Since the sampling of WMAP and Planck are different (see 
Sect. 3.2), we have used only the WMAP flux density measure- 
ments in the calculations here to avoid systematic errors. 

For many sources, especially at higher frequencies, we do 
not have good quality (i.e., 5 , > 3cr,) measurements for all years. 
We require at least four good data points (i.e., Nt„t > 4) at a given 
frequency to include a timeline in our variability analysis, and 
use S i - 3o-, as fiducial data points to fill in the years without 
good measurements. As discussed in Appendix A. 2, this leads 
to true lower Umits for;^'^ and Vrms- This approach also ensures 
that we have no variation in the sampling frequency of our data 
and formally have = 7 in all cases. We include all the avail- 
able Vnns values in Table 2 (expressed as a percentage), and flag 
strongly variable sources with a "V" and moderately variable 
sources with a "v" in the notes column. 

Table 3 summarizes, for each band, the total number of 
sources qualified for variability analysis, the number (and frac- 
tion) of strongly variable sources with their median fractional 
variability index Vmis (expressed as a percentage), and the num- 
ber (and fraction) of all the variable sources (both strongly and 
moderately variable) with their median Vtm^. The percentage of 
variable sources decreases with increasing frequency, while the 
level of variability, indicated by median V^ms, increases with 
frequency. A similar trend has been reported by Franzen et al. 
(2009) who suggested that it arises because sources tend to have 
larger flux density uncertainties and smaller flux densities at 
higher frequencies. Therefore only the sources with the strongest 
fractional variability can still be found to be variable at high 
significance at these frequencies. Our definition of Vtm& makes 
this explicit: (Vims) ( (o')w / (5')w) when calculated on a set 
of sources limited by variability significance (i.e., with a fixed 
threshold on x^ and hence Xims)- Since the flux density errors 
tend to increase and flux densities tend to decrease with fre- 
quency in our sample, (Vmis) will inevitably increase with fre- 
quency even when no true correlation of physical variability and 
frequency is present. To further test our hypothesis, we selected 
a subset of sources in our sample that are limited by a minimum 
^1,115 instead of Xims- We define Vrms so that all the sources (across 
all the frequencies) with Vmis ^ V'rms would have Xrms > 1 . We 
find Vtms - 22.7% in our sample, which leaves us a subsample 
of 56, 50, 31, and 20 sources at Ka, Q, V, and W bands, re- 
spectively. The median variability indices for this subsample are 
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Table 3. Variable Sources in our sample (198 sources in total). 



Strongly variable" Strongly + Moderately variable* 



Band 


Number 

(A's.r > 4) 


Number 


Median V,ms 
[%] 


Number 


Median V^ms 
[%] 


Ka 


165 


82 (49.7%) 


26.9 


97 (58.8%) 


25.1 


Q 


157 


67 (42.7%) 


27.2 


84 (53.5%) 


25.4 


V 


106 


32 (30.2%) 


31.4 


42 (39.6%) 


29.6 


w 


51 


15 (29.4%) 


33.4 


20 (39.2%) 


33.0 



Notes. These are sources that are variable at > 99% confidence level. '** These are sources that are variable at > 95% confidence level, including 
both strongly variable and moderately variable sources as defined in Sect. 4.1. 
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Fig. 4. Dependence of fractional variability index Vrms on mean 
spectral index (a) in Ka, Q, V, and W bands for variable sources. 
BLLs are shown by filled circles, FSRQs by open circles, and all 
the other sources by plus signs. The dashed lines indicate the 
level of Vrnis = 22.7%. 



found to be comparable in all four bands: 30.5% (Ka); 33.0% 
(Q); 31.4% (V); and 33.0% (W). Although our sample is not 
statistically complete, this confirms qualitatively the proposed 
explanation of Franzen et al. (2009) that the trend of increasing 
variability with frequency may be caused entirely by system- 
atics in source selection, and suggests that on the time scales 
and the frequency range we are probing, there is no significant 
increase of physical variability with frequency. Since different 
physical processes may dominate variability at short and long 
time scales (see Sect. 5), our finding is not necessarily incon- 
sistent with the conflicting results found on shorter time scales 
(e.g., Angelakis et al. 2012). 

We also tested for possible dependences of variability index 
on source spectral index and total flux density. For all the vari- 
able sources, we calculated the two-point spectral indices (i.e.. 



a 



Ka 



J, , uj,^, Cq, and Q-y ) for the years when good quality data 
points were available. Fig. 4 plots the fractional variability in- 
dex of each source against the time-averaged spectral index in 
each band. We use or^" to approximate the spectral index at Ka- 
band, Q-^^ for Q-band, for V-band, and for W-band. The 
variable sources are clearly dominated by flat-spectrum sources. 
Unlike Bolton et al. (2006) and Franzen et al. (2009), we do not 
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Fig. 5. Dependence of fractional variability index V„ns on mean 
flux density (Sy) in Ka, Q, V, and W bands for variable sources. 
BLLs are shown by filled circles, FSRQs by open circles, and 
all the other sources by plus signs. The dashed lines indicate the 
level of Vrms- The apparent trend of increasing variability with 
decreasing flux density is likely an artefact arising from the ap- 
plication of a significance cut in the presence of large measure- 
ment errors. 



find that sources with rising spectrum are more variable. We have 
also plotted the variability indices of all variable sources against 
their time-averaged flux densities at each frequency (Fig. 5). The 
apparent tendency for fainter sources to be more variable is not 
a physical trend, and is caused by the fact that fainter sources 
have larger relative flux density errors, and thus have a lower 
Vims when applying a significance cut (e.g., Xims > 1)- As can be 
seen from Fig. 5, when only sources with a fractional variabil- 
ity above Vims are considered, no significant correlation is found 
between the amplitude of variability and flux density. 

As mentioned in Sect. 2, most of the sources in our sam- 
ple can be classified as either luminous broad-line FSRQs or 
continuum-dominated BLLs. In Fig. 6 we plot the distribution 
of fractional variability index for the sources that are classified 
as variable in each band. There is no evidence that one class is 
more variable than the other. 
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Fig. 6. Histograms of fractional variability index for variable 
sources; FSRQs are shown in white and BLLs in grey. 



4.2. Modification of source number counts due to variability? 

The number of sources with good quality {S > 3cr) data across 
all seven years of WMAP observations is greatest at the WMAP 
Ka-band (33 GHz), with ~ 168 + 12 sources in each year The 
number decreases with frequency because most sources have flat 
or falling spectra, and flux uncertainties are also higher at higher 
frequencies due to the sensitivities of the maps. We therefore 
study the influence of source variability on number counts only 
in the Ka-band. In Fig. 7, we plot the differential number counts 
of the whole sample at the Ka-band in each year, including the 
ERCSC observation period. The ERCSC data points displayed 
here are again the ERCSC 30 GHz flux densities extrapolated 
to the Ka-band central frequency. Although they vary from year 
to year, our source counts are largely consistent with the pre- 
diction of extragalactic radio source evolutionary models from 
deZottietal. (2005) and Tucci et al. (2011). The departure of 
the data from the models suggests that our sample is incomplete 
below 2 Jy. We thus model the source count distribution dN/dS 
in Ka band as a power law KiS/Sy'f in each year, with a flux 
density cut of 2 Jy, using least-squares. The values we obtain are 
presented in Table 4. We apply a completeness correction to ac- 
count for the Galactic cut of \b\ > 5° that we used in selecting 
our sample. The slopes are close to Euclidean. It is clear that al- 
though the number counts fluctuate from year to year, they are 
consistent within the uncertainty. We therefore conclude that in- 
dividual source variability does not have a significant effect on 
the source number counts, i.e., the counts are statistically sta- 
ble, and so the contamination by unresolved sources in the CMB 
power spectrum that are estimated from the counts will also be 
stable. 



4.3. Correlation Analysis 

As mentioned in Sect. 1 (see also Sect. 5) geometric variability 
or low-peaked flares are likely responsible for the variability at 
the time scales and frequencies probed in our study. Both phys- 
ical processes predict that the average radio-millimetre spec- 
tral shape should remain nearly unchanged during flux varia- 
tions. We thus expect flux density variations in different bands 




S[Jyl 



Fig. 7. Euclidean normalized differential number counts at 
WMAP Ka-band (33 GHz) in each year The solid curve shows 
the total number counts of extragalactic radio sources predicted 
by the de Zotti et al. (2005) evolution model, and the dashed line 
represents the total number counts from the Tucci et al. (201 1) 
model. 

Table 4. WMAP Ka-band source counts in the flux density range 
2-20 Jy. 



Year 






sr-'] 




WMAPl 


86 


41± 


12 


-2.7+0.2 


WMAP2 


78 


29± 


9 


-2.5±0.2 


WMAP3 


81 


25± 


8 


-2.4±0.2 


WMAP4 


79 


25± 


8 


-2.4+0.2 


WMAP5 


78 


37± 


12 


-2.7+0.3 


WMAP6 


73 


36± 


12 


-2.8±0.3 


WMAP7 


84 


36± 


14 


-2.7±0.3 


ERCSC" 


75 


43± 


16 


-2.9±0.3 



Notes. Based on ERCSC 30 GHz flux densities extrapolated to the 
centre frequency of WMAP Ka-band. 



to be directly correlated. To quantify this correlation, we calcu- 
late the Pearson product correlation coefficient, r, between bands 
for each source. We note that this method quantifies the instanta- 
neous correlation between flux density changes, but is less sen- 
sitive to correlated changes with time lags between bands. The 
existence of significant time lags would reduce the correlation 
coefficient in our analysis. However, it is our intention in this 
study to distinguish between long-term (> 3 yr), frequency inde- 
pendent patterns, and short-term (< 1 yr), high-peaked flares (see 
Sect. 5.1). We do not employ more sophisticated methods such 
as discrete correlation function analysis to quantify possible time 
lags, as this is not meaningful for a sample containing only seven 
data points in each timeline. 

We analyse the correlations between Ka-band and each of 
the other bands (Q, V, and W). We restrict our analysis to sources 
which have been classified as strongly variable in the Ka-band, 
and satisfy the requirement of A^3o- > 4 and N2cr = 7 in each 
band; we do not require that the source to be classified as vari- 
able in the Q, V, or W bands. We include all data points with 
S > 2cr here to ensure that we have no missing data points in 
the time lines. Unlike in Sect. 4.1, we cannot use upper limits 
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for flux densities here as this would not lead to any useful lim- 
its on r, and the more restrictive requirement N-icr - 7 would 
leave us with too few time lines to analyse. We regard the vari- 
ability between two bands as "strongly correlated" if r > 0.875, 
which corresponds to a probability of less than ~ 1 % that the null 
hypothesis of an uncorrected variability is true (see Appendix 
A.3). We regard a value r < 0.3 as indication of no correlation, 
i.e., a 50% probability to reject perfect correlation. Sources with 
0.3 < r < 0.875 are considered "moderately" correlated. 

Our criteria select 55, 39, and 21 sources in the Q, V, and 
W band, respectively, for correlation analysis. Of these sources, 
58%, 46%, and 28% show strong correlation with the Ka band, 
and 40%, 49%, and 57% show moderate correlation. The major- 
ity of our sources show correlated variability between Ka and Q 
band, which confirms a similar finding by Franzen et al. (2009) 
between 15 and 32 GHz. Although there may be a moderate 
trend towards less correlated variability at higher frequencies, 
this should not be overstated as it may be due to the larger error 
bars aff'ecting these bands. In general we find that very few vari- 
able sources show no direct correlation between Ka and higher- 
frequency bands. This result may suggest that variability in these 
frequencies on long time scales is less frequency-dependent than 
on short time scales (e.g., Angelakis et al. 2012). 

5. Variability patterns in blazars 

In this section we describe the variability patterns expected in 
simple models motivated by current understanding of AGNs and 
extragalactic jets, and search for the signatures of these patterns 
in the time series of sources in our sample. 

5. 1 . Physics of blazars and models for variability 

The term blazar was introduced in the late 1970s to denote 
the combined class of optical violent variable (OVV) quasars 
and BL Lac objects (Angel & Stockman 1980). In the radio-to- 
millimetre range, blazars show featureless, often highly polar- 
ized, continuum spectra which are attributed to partially self- 
absorbed synchrotron emission from relativistically boosted, 
compact jets emerging from AGN at a small angle 6 to the line 
of sight (e.g., Begelman et al. 1984). If T = 1/ ^J]~^ » 1 
is the bulk Lorentz factor of the relativistic jet plasma, the 
Doppler boosting factor D - [r(l -jScosfl)]"' is ~r for 
6 < F This causes the flux density of the synchrotron emis- 
sion to dominate over unboosted thermal components of the 
AGN, as Sv oc D^^". The Doppler factor, bulk Lorentz fac- 
tor and inclination angle of the jets can be estimated by sev- 
eral methods. Firstly, by comparing number counts of blazars 
with those of unboosted radio galaxies (Urry & Padovani 1995). 
Secondly, from superluminal apparent speeds of VLBl compo- 
nents (e.g., Blandford et al. 1977; Ghisellini et al. 1993; Zensus 
1997; Kellermann et al. 2004). And thirdly, by exploiting physi- 
cal constraints on models of radiative emission (Ghisellini et al. 
1998; Lahteenmaki et al. 1999; Lahteenmaki & Valtaoja 1999; 
Hovatta et al. 2009). All these methods allow a consistent esti- 
mate of D ~ r ~ 10 and 9 ~ 5°. Estimates of individual Doppler 
factors and viewing angles, however, are too model-dependent 
to be considered reliable (e.g., Hovatta et al. 2009). 

Blazars are strongly variable down to time scales of ~ 1 day 
or less. Such rapid variability is usually attributed to the pres- 
ence of shocks in the jets, which are expected to arise from in- 
stabilities inside the jet or from interactions of the jet with an ex- 
ternal medium (Begelman et al. 1984). Several phenomenologi- 



cal approaches to the formation and evolution of shocks in jets 
have been suggested, ranging from persistent moving or stand- 
ing shock waves (Marscher & Gear 1985; Marscher et al. 2008) 
to rapid sequences of shocks formed over a large range of scales 
in the jet (Spadaetal. 2001; Guettaetal. 2004; Rachen et al. 
2010). Shocks in AGN jets can also provide the sites for effi- 
cient particle acceleration (e.g., Biermann & Strittmatter 1987) 
needed to produce the observed high energy gamma-ray emis- 
sion (see, e.g., Bottcher 2012). 

Another type of variability, usually called geometric 
variability, is expected in a number of phenomenolog- 
ical models describing helical or precessing motion in 
jets (Camenzind & Krockenberger 1992; Steffen et al. 1995; 
Valtaoja et al. 1989; Wagner et al. 1995; Viflata & Raiteri 1999; 
Ostorero et al. 2004; Valtonen et al. 2009, 201 1). Because of the 
strong dependence of the beamed flux on the Doppler factor 
(oc D^^"), small changes of the inclination angle or the bulk 
Lorentz factor F of the jet will induce large changes in the ob- 
served flux. For example, in a source with F = 10, 6* = 5°, and 
Q' = 1 , a change of inclination angle by one degree towards the 
observer increases the measured flux by a factor of two. 

These two types of variations are expected to affect the 
spectrum of the source diff'erently: while shocks naturally pro- 
duce spectral changes by injecting a new population of ener- 
getic particles, geometric effects should not have a large effect 
on the emitted spectrum. Detailed modelling, however, shows 
that this simple criterion has to be used with care. For example, 
the strongly inverted "edges" often seen in the radio-millimetre 
spectra of variable blazars (Valtaoja et al. 1988; Brown et al. 
1989; Krichbaum et al. 2008; Planck CoUaboration XIV 2011; 
Planck Collaboration XV 2011) can be explained both by syn- 
chrotron self-absorbed emission from a compact, shocked region 
with temporarily enhanced emission compared to the surround- 
ing jet (e.g., Valtaoja et al. 1988; Brown et al. 1989; Tiirler et al. 
2000; Rachen et al. 2010), and by decreasing inclination an- 
gle in a steady but rotating helical jet (Ostorero et al. 2004). 
The same ambiguity affects the interpretation of helically mov- 
ing VLBl components (Savolainen et al. 2002; Bach et al. 2006; 
Rastorgueva et al. 2009; Jorstad et al. 2010). On the other hand, 
achromatic variability, canonically expected from geometric ef- 
fects, can also be explained by shock-induced flares if they are 
observed only above their peak-frequency (so-called low-peaked 
flares, Valtaoja et al. 1988). Thus, although the overall proper- 
ties of variability from shocks and from geometric effects can be 
quite different, the spectral behavior during flares may not be a 
rehable way to distinguish the processes. 

5.2. Short-term vs. long-term variability 

A major difference between geometric and shock induced vari- 
ability is the observed time scale. In a relativistically boosted 
emitter, all intrinsic time scales of variability are compressed by 
the Doppler factor, Tobs - TintlD <K Tim. This applies in particu- 
lar to the time scales connected to the formation and propagation 
of shocks in the jet, and allows emission regions of sub-parsec 
size to produce variability on observed time scales down to a few 
days. In contrast, the physical time scale of geometric changes 
are not Doppler compressed, as the rotation is transverse to the 
jet propagation (robs ~ Tint). If indeed both processes contribute 
to variability with comparable intrinsic time scale and ampli- 
tude, geometric variability should then dominate the amplitude 
on long time scales. We illustrate below how this affects the pat- 
terns observed in long-term averaged data. 
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Let Q{t) and T(t, v) denote the time structure of geometric 
and shock induced (flare) variability respectively, with the total 
flux of the source being S(_t, v) - @{t) T(,t, v). The WMAP data 
used in this paper provide regularly sampled 1-year averages, 



I.e., 



' WMAP' 



(ti, Vb) 



dtg(t) T{t,vb) w(ti-t), 



(8) 



where f, is the end-time of operation year ;, v;, is the central fre- 
quency of the appropriate WMAP band, and w is the exposure 
function of WMAP in this band (which we assume is approx- 
imately the same in every operation year). Eq. 8 describes a 
low-pass filter, suppressing all fluctuations in S (t, Vb) on time 
scales T <sc Ar = f; - f,„i ~ 1 yr. This is typically the case 
for the flare process !F(f, v^,), but not for the geometric process 
Q{t), as its characteristic time scale is not Doppler compressed. 
Therefore, although large amplitudes of 'Fit, Vh) may cause an 
irregular structure in S (f, v/,), characteristic patterns of ^(f) may 
become more evident when using averaged data. We note that 
these considerations apply to the WMAP data, but not to the 
Planck ERCSC data in general. Most ERCSC flux densities rep- 
resent averages of one or two snapshots, except for a few sources 
near the ecliptic poles which were observed more often. The flux 
densities thus retain a significant fraction of the short-term vari- 
ability amplitude. We therefore discuss the WMAP and ERCSC 
data separately. 

5.3. Test for variability patterns 

Geometric variability is not bound to specific simple patterns, as 
the direction of jet emission can change in potentially chaotic 
ways. However, a major motivation for proposing geometric 
variability is the possible connection to helical jet structures 
shown by VLBI monitoring (see Sect. 7.1), which are thought 
to originate in a rotating jet base. As discussed in Appendix C, 
a simple rotating relativistic jet will produce a flux variability 
pattern of the form 



5rot(0 = 



(1 -I- Acos[w(f- fo)]) 



3+a 



(9) 



Assuming the amplitude A and frequency at of the oscillation are 
time-independent, this model involves four free parameters (the 
spectral index a is determined from the data). 

A similarly simple phenomenological description of the time 
evolution of flares caused by shocks has been introduced by 
Valtaoja et al. (1999, hereafter VLTL), based on observed flare 
properties and a generalization of the "shock in jet" model 
(Marscher & Gear 1985; Valtaoja et al. 1988). It has been suc- 
cessfully applied to decompose single-frequency blazar light- 
curves into a constant background and a number of separate 
exponential flares (see also Nieppola et al. 2009). In the orig- 
inal VLTL model, a flare is described by an exponential rise 
with time constant t, followed by an exponential decay with 
time constant Tf ^ 1.3ti-. The time scales of such flares can 
vary between a few months and ~ lOyr, with a median of ~ 3 yr 
(Nieppola et al. 2009), which is comparable to the time scales 
probed in this paper This motivates a four-parameter VLTL pat- 
tern as a heuristic description of a single flare or activity phase: 



ma> 



.('m»-f)/l-3T 



t<tn 
t>t„ 



(10) 



Both models are clearly simplified. Realistic models of 
helical jets involve complex dependencies on more param- 
eters (Camenzind & Krockenberger 1992; Steffen et al. 1995; 



Villata & Raiteri 1999), and may potentially include static in- 
homgenities in the jets (Ostorero et al. 2004). For the VLTL 
model, Nieppola et al. (2009) showed that the relation Tf ^ LSt^ 
holds only as a gross average, while the individual flare decay- 
to-rise time scale ratio varies between 0.3 and 5. Moreover, it is 
likely that more than one flare appears during the decade probed 
by our data. We therefore emphasize that our quest is not to pro- 
vide fits to the data with realistic models, but to test how well 
our data can be represented by the simplified models discussed 
here. 

We limit our analysis to sources that show strong corre- 
lated variability between Ka band and at least one other band 
(Sect. 4.3), which ensures a better detection of common patterns 
between bands. In order to use information from all the bands to 
determine the best fit parameters of a pattern and its fit quality, 
we normaUze all the data points to a common level at 5 - S 
where 



r 



2(5^) 



(11) 



This normalization is somewhat arbitrary, and was chosen 
mainly for comparability with external data (see Sect. 5.4.3). 
Using an initial guess for both the VLTL and rotating patterns, 
we minimize the residual with respect to the model using the 
Nelder-Meade simplex method implemented in the GNU scien- 
tific library^. We then calculate the residual Xres, defined as 



IN- 



(12) 



where p = 4 is the number of fitted parameters. We accept a 
model fit if X,es < 1, or if it reduces the mean of Xn„s (defined 
as in Eq. 4) from all the bands by a significant factor, i.e., X^es < 

3 ^rms ■ 

The one-year averaged WMAP flux densities have equal 
weighting in each year, and therefore it is correct to space the 
points by exactly one year when fitting models, even if the exact 
epoch of each (and hence the reference time of the fit) is uncer- 
tain by a few months. We do not use the Planck ERCSC flux 
densities in the fitting as they have different weighting, and the 
time spacing relative to the WMAP observations is difficult to 
determine. 



5.4. Results 

Among the 32 sources in our sample that show strong corre- 
lated variability between Ka-band and at least another band, we 
find 19 (60%) whose long-term averaged variability pattern is 
well represented by a simple four-parameter model. The other 
13 sources show more complex variability patterns. 



5.4.1. Simple Patterns 

A preference for a rotating jet pattern is found for 
sources J0428-3756, J0457-2324, J0721+7120, Jl 159+2914, 
J1310+3220, J1337-1257, J1733-1305, J1924-2914, and 
J2235-4835. Most of these sources show strongly correlated 
variability in more than two bands. The rotation periods found 
are in the range 3-5 yr The shortest rotation period of 3 yr is 
found for J1337-1257, which showed more than two complete 
oscillations within our observing period. The WMAP/Planck 
light curves of these sources, as well as the fitted pattern, are 



http://www.gnu.org/software/gsl/ 
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Fig. 8. Upper panels: WMAP/Planck light curves for the sources that are better fit by the rotating jet model. The WMAP data are 
shown by filled symbols and the later Planck data are shown by open symbols. Lower panels: Residuals between the normalized 
WMAP flux densities (Eq. 11) and the rotating jet model. Grey dash-dot lines indicate the zero level. 



shown in Fig. 8. Although the fitting is done on the scaled 
WMAP flux densities, we plot the original WMAP flux densi- 
ties to demonstrate that the model traces the variation trend of 
the original data very well. The residuals shown in the bottom 
panels are, however, the differences between the normalized flux 
densities and the fit. 

For sources J0106-4034, J0253-5441, J0403-3605, 
J0423-0120, J0440-4333, J0530+1332, J1613+3412, 
J1635+3807, J1743-0350, and J2258-2757, both the ro- 
tating jet model and the VLTL-flare model yield a reasonable 
match to the long-term flux variation pattern. Fig. 9 gives their 



WMAP/Planck light curves, together with the pattern fit. With 
the exception of J0423-0120, the rotation periods of these 
sources in the rotating jet model are found to be > 7 yr, and the 
flare time scales in the VLTL-flare model are ~ 0.5-2 yr For 
five sources with rotation periods > 10 yr (r > 1 .4 yr) we mostly 
see only the increasing or decreasing part of the pattern in the 
time window of our observations. 

The Planck ERCSC flux densities, though not used in the fit- 
ting, are found to be in good agreement with the model fits in 
some cases, e.g., J0457-2324, J1310+3220. When both models 
fit the WMAP data, the Planck flux densities sometimes prefer 
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Fig. 9. Top panels: WMAP/Planck light curves for the sources that are well fit by both the VLTL-flare model and the rotating jet 
model. The WMAP data are shown by filled symbols and the later Planck data are shown by open symbols. Middle panels: Residuals 
between the normalized WMAP flux densities (Eq. 11) and the VLTL-flare model. Bottom panels: Residuals between the normalized 
WMAP flux densities and the rotating jet model. Grey dash dot lines indicate the zero level in the residual plots. 
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Fig. 9. Continued. 



one model over the other, e.g., J0403-3605 for the rotating-jet 
pattern, and J2258-2757 for the VLTL-flare pattern. In about 
50% of the cases, the Planck data are in disagreement with the 
fitted pattern. This is likely due to the unmodelled short-term 
variability and the difference between the actual Planck obser- 
vation times and the referenced time (set to be exactly 2 yr later 
than the last WMAP data point). 

5.4.2. Sources with complex variability 

There are 13 sources for which neither the VLTL-flare model 
nor the rotating jet model fits the variation pattern well. These 
sources can be divided into two classes: (1) those showing 
highly-correlated variability across all bands; and (2) those 
showing significant diff'erences between the variability patterns 
at diff'erent bands. 

In the first category, we find bright FSRQs like 3C273, 
3C279, and J0538-4405 (Fig. 10). These sources would proba- 
bly be fitted by a correlated pattern comprising several VLTL 
flares, or a rotating pattern plus long-term flares. Sources in 
the second category all show moderately correlated, partially 
inverted spectra at higher frequencies. A famous example is 
3C 454.3, for which our data capture the well-studied millimetre 
flare from 2005 (Krichbaum et al. 2008). A very strong inversion 
in the V-W spectral index is also found in source J2232-(-1143 
between 2002 and 2004. In addition, the light curves of BL Lac, 
3C446, J1058-H0134, and JI549H-0236 all display similar spec- 
tral behaviour (see Fig. 11). Such signatures of spectral depen- 
dence of variability are well known at shorter time scales, and 
are expected from physical models of both shock-induced flares 
and geometric effects in rotating helical jets (Villata & Raiteri 
1999; Ostorero et al. 2004). 

5.4.3. Comparison with long-term monitoring data 

A limitation of our analysis is the small dynamic range in time 
scales within which our method can distinguish between differ- 
ent patterns; time scales smaller than one year are not resolved, 
and for time scales longer than a few years our fits become am- 
biguous, as discussed above. 



For the three sources J0423-0120, J0721+7120, and 
Jl 3 10-1-3220 that are included in the 37 GHz long-term monitor- 
ing program on the Metsahovi radio telescope (A. Lahteenmaki, 
private communication, see also Teraesranta et al. 1998), we can 
test the validity of our finding over a long time period (> 20 yr). 
In order to make the Metsahovi data comparable to our WMAP 
data, we averaged them into one year bins (see Appendix B). 
Fig. 12 shows the the WMAP Ka and Q band light curves along 
with the binned Metsahovi 37 GHz data. In all three cases, we 
find excellent agreement between the Metsahovi binned data and 
the WMAP data. This confirms that long-term averaged data pro- 
vide a good and consistent representation of the long-term vari- 
ability trends, regardless of the details of the averaging method. 

The three examples both highlight the potential strengths 
and reveal the various problems of our approach. For source 
J0423-0120, neither the Metsahovi data nor the Planck data sup- 
port the rotating-jet pattern that is compatible with the WMAP 
data; rather, the binned Metsahovi data suggest that the source 
underwent a single, strong activity phase between 2001 and 
2005, weU fitted by a VLTL flare profile with t ^ 0.5 yr, 
superimposed on a flux density which remained roughly con- 
stant between 1980 and 2010. For source J0721-I-7120, both the 
Metsahovi and the Planck data are consistent with the fitted rota- 
tional pattern between 2002 and 2010, but before 1998 the long- 
term average flux density of the source was continuously in a low 
state. This could be interpreted as the onset of a new rotating 
pattern as expected, e.g., in lighthouse scenarios (Steffen et al. 
1995). An alternative explanation would be to assume two con- 
secutive long-term flares (or flaring phases) of r ~ 1 yr in the 
years after 2002. Finally, for source J1310+3220, the Metsahovi 
and Planck data are consistent with a continuation of the rotat- 
ing pattern fitted to the WMAP data alone at least for the period 
1998-2010, covering almost two full oscillations. Additionally, 
the long-term averaged 37 GHz light curve observed between 
1983 and 1990 fits the same pattern, with a consistent phase, 
as expected of geometric variability from stable rotating jet; the 
departure from the pattern between 1990 and 1998 could be ex- 
plained by a single superimposed flare (t ~ 1 yr). 

6. Variability in unbeamed sources 

Our sample includes some bright radio galaxies, which in the 
standard unification model (e.g., Urry & Padovani 1995) are the 
unbeamed counterparts of blazars. We note that for a source 
with inclination angle » F the contribution of the compact, 
highly relativistic jets is not only unbeamed, but deboosted with 
a Doppler factor D ~ 0.1 if the source is exactly viewed from the 
side. This suppresses its contribution to the total emission of the 
source by a factor > 10* relative to a beamed blazar. Therefore, 
in such sources all the emission must come from the extended 
jets and lobes, which are subdominant in highly beamed blazars. 

The extended radio galaxies in our sample, most notably 
Cen A, Cyg A, and M 87, are all classified as non-variable at all 
frequencies, which is expected because the emission regions in 
these sources extend over several kiloparsecs, too large to show 
significant variability on the time scales we probe. This is a con- 
sistency check on the reliability of our photometry. 

6.1. A millimetre-wavelength fiare in ttie giant radio galaxy 
Pictor A 

The giant radio galaxy Pic A, however, is an exception: it shows 
moderate to strong variability in the Ka, Q, and V bands, with 
values of 17.0, 18.3, and 15.6, respectively. The flux densities 
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Fig. 10. WMAP/Planck light curves for the sources J0538-4405, 3C 273 and 3C 279. The WMAP flux densities are shown by filled 
symbols and the Planck flux densities are shown by open symbols of the same shape. The dashed lines indicate the 3cr level at each 
band if in the plotted range. 
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Fig. 11. WMAP/Planck light curves for the sources J1058+0134, J1549+0236, BLLac, 3C446, J2232+1143 and 3C 454.3. The 
WMAP flux densities are shown by filled symbols and the Planck flux densities are shown by open symbols of the same shape. The 
dashed lines indicate the 3cr level at each band if in the plotted range. 



in the WMAP Ka, Q, and V band were generally increasing until 
2008. Planck ERCSC data suggest that in 2010 the flux density 
returned to its 2002 value for all bands (Fig. 13). Thex^ values 
suggest a combined probability of ~ 10"^ for constant flux den- 
sity in all bands, so the variation is highly significant. From the 
shape of the variation, which is essentially coherent in all bands, 
we can estimate an observed flare time scale of ~ 3 yr, assuming 
an exponential rise, and a somewhat shorter decay time scale. 

Pic A is a Fanaroff-Riley class II (FR II) radio galaxy at a red- 
shift z - 0.035. It shows two prominent radio lobes, with bright 
hot spots terminating the jets at distances of about 140 kpc from 
the central AGN. They are connected to the AGN by extended 
jets, which appear to propagate with mildly relativistic speeds 
(Meisenheimeret al. 1989). Within these central jets, at a dis- 
tance of ~ 30 kpc from the AGN, Chandra observations provided 
evidence for a transient X-ray knot, which apparently decayed 
within the time scale of ~ 1 yr (Marshall et al. 2010). Assuming 
that the X-ray radiation was produced by synchrotron radiat- 



ing electrons, Marshall et al. (2010) concluded that the magnetic 
field in the emission region must have been ~ 2 mG, about 100 
times larger than the canonical estimates for the jets of radio 
galaxies. This magnetic field, however, is still too weak to ex- 
plain the time scales of the flare observed by WMAP at mil- 
limetre wavelengths (also assuming synchrotron origin), there- 
fore we consider it unlikely that the two flares originate in the 
same region of the jet. 

An alternative explanation is that the increase in millimetre 
flux density is not from the extended jet, but from the central flat- 
spectrum core, which is almost as bright as the dominant west- 
ern hot spot at 20 GHz (Perley et al. 1997; Burke-Spolaor et al. 
2009; Sajina et al. 201 1). Interpreting the variable component as 
unboosted (D ~ I) emission from a central blazar, and assuming 
standard jet parameters, the emission region would have an in- 
clination angle to the line of sight of 9 ^ 25°, and the observed 
factor-of-two increase in flux density within three years could 
be explained by a change of inclination angle by » 3° towards 
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Fig. 12. WMAP/Planck light curves for the sources J0423-0120, J0721+7120, and J1310+3220 at 33 GHz {black filled circles) 
and 41 GHz (red filled triangles), with the Metsahovi 37 GHz long-term monitoring data added {grey stars). For J1310-(-3220, the 
residuals between the Metsahovi data and the best fit model are shown. 
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Fig. 13. WMAP/Planck light curves of the FRII radio galaxy 
Pic A. 



the observer. However, this would require a rotational speed sig- 
nificantly larger than that inferred from VLBI observations of 
the central jet of Pic A, if these are interpreted in a precessing 



jet model (Tingay et al. 2000). An interpretation of the flare in 
Pic A as due to a shock in the jet would be more consistent with 
the msialigned-blazar picture: with a Doppler factor D ~ 10 the 
Pic A flare would be comparable in intrinsic time scale and peak 
frequency with the exceptional outburst in 3C 454.3 observed in 
2005. 



7. Discussion 

7.1. Geometric variability and lielical jets 

Rotating jet patterns have been discussed in the literature in con- 
nection with VLBI observations of morphologically "curved" 
or "bent" jets, and of apparent helical motion of parsec-scale 
jet components (Zensus 1997; Kellermann et al. 2004). Among 
the sources for which we find a variability pattern matched 
with the rotating jet model, curved and potentially helical 
VLBI jets have been observed in J0423-0120 (Wagner et al. 
1995; Britzenetal. 2001), J0530+1332 (Pohl et al. 1996; 
Caietal. 2006), J0721H-7120 (Bach et al. 2005), J1159-H2914 
(Hongetal. 2004; Zhao et al. 2011), J1337-1257 (Lister etal. 
2009), J1613-H3412 (Piner & Kingham 1997), and J1924-2914 
(Shen et al. 1999; Lister et al. 2009). Helical structures have also 
been found in the VLBI jet of FSRQ J13 104-3220 (Cassaro et al. 
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2002), for which a rotating pattern with essentially constant pe- 
riod is largely consistent with the flux variation pattern shown in 
the Metsahovi 37 GHz data (Sect. 5.4.3). 

Although our data are consistent with geometric variabil- 
ity, all the above-mentioned sources also show flaring behaviour 
on short time scales. An extreme example is the BLLac object 
J0721-)-7120, which showed rapid variability on time scales of a 
day or less (Ostorero et al. 2006). In some sources this flaring ac- 
tivity may have obscured an underlying rotating jet pattern, or in 
some cases our geometric model may be too simple. For exam- 
ple, BL Lac shows complex long-term variability in our data, but 
its variations have been sucessfully modelled over much longer 
time scales as a purely geometric effect, assuming a rotating in- 
homogeneous helical jet (Bach et al. 2006; Villata et al. 2009). 

Long-term rotation cycles have been identified in 
VLB! observations of 3C 454.3 (Qian et al. 2007) and 
3C279 (Carrara et al. 1993); both sources show strong variabil- 
ity correlation between bands, but a complex variability pattern 
in our data (Sect. 5.4.2). Source J0423-0120 shows a pattern 
compatible with the rotating jet model in the WMAP data, but 
the 25-year Metsahovi 37 GHz light curve does not support a 
rotational pattern. 

There are a few helical jet candidates that are not included 
in our variability pattern test due to weak correlation of vari- 
ability between WMAP bands. These include 3C 84 (NGC 1275, 
Krichbaum et al. 1992), 3C345 (Steff^en et al. 1995; Qian et al. 
2009), J 1800+7828 (Cassaro et al. 2002), and the famous binary 
black hole candidate OJ 287 (Siflanpaa et al. 1988; Villata et al. 
1998; Valtonenet al. 2009, 2011). The FSRQ J1833-2103 is a 
gravitationally lensed rotating jet source (Guirado et al. 1999; 
Nair et al. 2005) that is also excluded from our pattern test, ow- 
ing to the lack of correlation between flux variations in the W- 
band and the other WMAP bands, although visual inspection sug- 
gests a rotating -jet pattern with a period of < 10 yr 

In summary, rotating jets are common features in AGN, but 
only some of them display a characteristic sinusoidal pattern in 
the long-term averaged WMAP light curves. 

7.2. Relation to other frequency bands 

The long-term variability and possible periodic behaviour of 
AGN have been studied using many methods (e.g., structure 
function, wavelet, periodogram, DCCF; see Sect. 5.1). Evidence 
for periodic behaviour has mostly come from optical data, in 
particular for bright sources where archival optical data provide 
time series extending over many decades (Villata et al. 1998; 
Fan & Lin 2000). These periodicities have generally been inter- 
preted as resulting from geometric effects, potentially in con- 
nection with shocks (Lainela et al. 1999). Whilst large samples 
of radio time series spanning more than 20 yr have not provided 
strong support for periodicities, WEBT campaigns on BL Lac 
and AO0235+16 (J0238+1637) have shown that time lags be- 
tween radio and optical variability can be explained as a geomet- 
ric effect, where rotation successively moves regions emitting at 
different frequencies into positions of maximal Doppler boosting 
(Ostorero et al. 2004; Bach et al. 2006; Villata et al. 2009). 

Gamma-ray emission in blazars is likely to a have a differ- 
ent origin from the radio emission (e.g., Ghisellini et al. 1998; 
Rachen 2000; Bottcher 2012), so comparison of the long-term 
variability in radio and gamma-ray data could provide an impor- 
tant cross-check for the geometric interpretation, which predicts 
that these data should show long-term correlations like those 
seen in the radio and optical data. Among our sources, 133 are 
present in the second source catalogue from the Fermi Large 



Area Telescope (2FGL, Ackermann et al. 2011), 105 of which 
show significant variability in gamma-rays. The uniform sam- 
pling of Fermi-GST makes it an ideal instrument for such study. 

7.3. Caveats 

We find simple, sinusoidal patterns on time scales < 10 yr in only 
nine out of 198 sources, i.e., 4.5% of the full sample. This is con- 
sistent with the results of Ciaramella et al. (2004), who found pe- 
riodic behaviour in about 4% of their irregularly sampled blazar 
light curves spanning several decades. For most sources, the 
time scales are too long or the variability is too complex to fit 
our simple models. More complex, physically motivated mod- 
els are possible, and have been used to model selected blazar 
light curves. However, the added flexibility of the more complex 
models allows both geometric models (Villata & Raiteri 1999; 
Ostorero et al. 2004; Bach et al. 2006) and shock-based mod- 
els (Hovatta et al. 2008; Nieppola et al. 2009) to fit the data. 
For example, the exceptional outburst observed during 2005 
in 3C 454.3, has been modelled successfully in both geomet- 
ric (Villata et al. 2007) and shock-based scenarios (Rachen et al. 
2010). One approach to discriminate between models of dif- 
ferent complexity is to adopt Bayesian classification methods 
(Rachen 2013). It is possible that such methods, in conjunction 
with more data and better models, may eventually elucidate the 
true cause of long-term variability in blazars. 

The present data merely demonstrate that a surprisingly sim- 
ple sinusoidal pattern, indicative of geometric variability, is a 
good representation of the long-term averaged variability be- 
haviour for more than half of the sources in our selected sam- 
ple of 32 sources. It is not yet possible to exclude any of the 
more complex explanations for blazar variability that have been 
proposed. 

8. Conclusions 

We have used data obtained from the WMAP and Planck satel- 
lites to analyse the long-term flux-density variations of a sample 
of extragalactic radio sources selected from the Planck ERCSC. 
We have used flux densities averaged over one-year periods, 
which suppresses short-term variability and makes the data set 
particularly suitable for studying underlying, long-term variabil- 
ity patterns. 

At WMAP Ka, Q, V, and W band, we found 82, 67, 32, and 
15 sources to be variable at > 99% confidence level. The number 
of variable sources decreases at higher frequencies owing to the 
higher flux density uncertainties at these frequencies. The am- 
plitudes of variation are found to be comparable in the different 
bands, and are not correlated with either the flux densities or the 
spectral indices of the sources. At Ka band, we modelled the 
source number counts in each year, and found that although the 
fitted parameters for the model vary from year to year as a con- 
sequence of individual source variability, the variation is within 
the uncertainty of the fitted parameters, and therefore the num- 
ber counts derived from the source sample remain stable over 
the years. Almost all of our sources show correlated variability 
between bands. Essentially all of the sources in our sample are 
radio loud AGNs, and most of them are blazars. We do not find 
any significant difference between the two subclasses of blazars, 
FSRQs and BLLs. 

As expected from theory and VLBI observations, both shock 
induced flares and geometric variability due to helical or pre- 
cession jet motion may contribute to the total variability at long 
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time scales. In more than half of the sources that show strong 
correlated variability between Ka and another band, the long- 
term averaged variability shows the sinusoidal patterns expected 
from simple jet rotation, with periods greater than a; 3 yr. A sim- 
ple exponential flare model is less successful in representing our 
data, but definite conclusions on the nature of variability cannot 
be drawn because more complex models cannot be ruled out. 

Most extended radio galaxies in our sample show no sign of 
variability. An exception is the extended radio galaxy Pic A, for 
which we found significant changes in the flux density between 
2002 and 2010, which could arise from the central "misaligned 
blazar" nucleus if the inner jet has a small angle to the line of 
sight. 

The unexpected observation of a millimetre-flare in an ex- 
tended radio galaxy underlines the potential of cosmology 
probes and other multi-survey instruments for research on vari- 
ability and transient phenomena, as they allow comparisons be- 
tween many, identically sampled full sky surveys. 
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Appendix A: Quantifying variability in astronomical 
data sets 

A.1. Statistical significance and variability index 

To quantify significant changes in the measured flux density of 
a source, recent astronomical literature has adopted a so-called 
de-biased variability index. 



<5) V A' 



(A.l) 



The expression is usually attributed to Akritas & Bershady 
(1996), but the form of Eq. A.l was first applied to time series 
by Sadler et al. (2006). For the special case that all measurement 
errors are equal, cr,- - cr, Vi-^s can be written as 



(A.2) 



Assuming is large and the bias through diminishing one de- 
gree of freedom is small, this equation relates to the statistical 
significance of variability in a plausible way: it shows that for 
X^/N > 1, i.e., an average scatter in excess of the measurements 
error. Vims is positive, and is given as a fraction of the average 
total flux density by scaling with cr/{S). Fot x^/N < 1, Vim,s 
becomes imaginary, indicating that the data scatter is less than 
expected. Data sets showing a large fraction of sources with 
imaginary V^ms are likely to have overestimated the measure- 
ment errors and thus underestimated the fraction of truly variable 
sources. 

From the derivation of this expression in Akritas & Bershady 
(1996), we note that, in spite of the notation of individual er- 
rors, Eq. A. 1 is based on a valid statistical estimator only for 



the case of equal errors. Moreover, it has been derived under an 
assumption of large sample sizes (see also Franzen et al. 2009). 
Comparing with Eq. A.2, we immediately see where this con- 
straint comes from: (a) the failure to reduce the degrees of free- 
dom in the normalized and (b) the use of the arithmetic mean 
(t/{S), which does not minimize except in case of equal er- 
rors. Correcting for this, i.e., doing the transformations 



±1 

N 



X" 



(A.3) 



and introducing {(t)^ as an estimate for the average error in the 
scahng factor, we obtain Eq. 4. 

A.2. Dealing with incomplete data 

In this paper, we are limited not only by the small sample size, 
but also by the small fraction of above 3>cr detections in a time 
line in many cases. However, given that enough "good" (> 3cr) 
data are available to achieve a reasonable estimate for the mean 
source flux density, and upper limits in the years in between can 
provide valuable information on the variability of the source, we 
chose to treat all upper limits as data points with S ,■ = 3crnoisg and 
o"/ - o"noise, and use Eq. 7 to calculate Since a reasonable es- 
timate on the average source flux density S above the noise level 
is vital for a meaningful analysis, we require more good data 
than upper limits in order to include a time line in our analysis, 
i.e., N-icr > 4. We also note that the value obtained in this way 
for variability strength and significance is a lower limit, as the 
distribution is cut on one end by our noise level, which reduces 
the scatter 

A.3. Pearson product correlation coefficient 

For our correlation analysis, we calculated the Pearson product 
coiTelation coefficient, defined as 



r(Ka,X) 



1 



VvarCSK^ )\ Vvar(5^) 



(A.4) 



between the Ka band and other bands X = Q, V, W showing 
significant variability. As all errors in the same band for the same 
source are approximately equal, we took {S ) and var(5 ) as the 
straight, unweighted arithmetic mean and variance of the data in 
each time line. For approximate Gaussian distributions of S in 
each band, and large enough samples, the variable 



t - r 



N-2 



follows Student's f-distribution for A^-2 degrees of freedom. We 
note that none of the above assumptions really holds for our data 
sets, but more exact derivations of correlation confidence do not 
help as they also depend on the assumption of Gaussianity. 

Appendix B: Binning and de-biasing Metsahovi 
monitoring data 

To assure compatibility with the WMAP results, we need to bin 
the Metsahovi data to match WMAP operation years, starting 
from August of a given year through July of the next year We 
have to consider that, unlike WMAP with its fixed scanning pat- 
tern, monitoring programs like Metsahovi do not sample data in 
an unbiased way. Telescope time is usually focused on sources 
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that show rapid flaring; sources that are not in an active phase are 
observed less frequently. Therefore, we have more data points in 
high states than in low states of a source, leading to strong bias 
in simple averages. To reduce the bias, we bin the data using a 
weighted average 

where j denotes the number of the bin, and the sum runs over all 
data points included in the bin with individual flux densities S p. 
Here, Af, - r,+i - f, represents the time diff'erence between ad- 
jacent flux density measurements in the monitoring, and is used 
as the weight factor in the average to compensate the bias in- 
troduced by the irregular sampling rate. In a similar way, we 
determine a weighted average of squared flux density {S^}j, and 
obtain the standard deviation 

= ^|(S~)J-{S}], (B.2) 

which describes the amount of variation of the source flux den- 
sity within the one-year bin. We note that this should not be con- 
fused with a measurement error; rather it reflects the standard 
deviation of the distribution of data points within the bin. It may 
therefore be regarded as a measure of the strength of short-term 
variability during the respective one-year period. 



Appendix C: Brightness variation in relativistically 
beamed rotating emitters 

We consider a beamed emitter rotating with angular velocity co 
around a cone of opening angle O about a central axis of a jet. 
Let be the angle of the line of sight to the central axis of the 
cone. The inclination angle of the jet at time t to the line of sight 
is then given by 

cos 9 = cos <J> cos + sin <1> sin cos ojt , (C.l) 

which can be simply written as a+b cos ojt. For a beamed emitter, 
fl » 1 and b <s: 1. Inserting this into the definition of the Doppler 
factor and setting a - 1 , we obtain 



1 + A cos cot 

with A ^ b/3/(l - (3) < I, where we used 1 -yS^ = l/T^, and 
eventually set/? = 1. As 5,, D^*^", this leads to Eq. 9. 

If we consider a combined variability process which com- 
prises both flare and geometric variability (see Sect. 5.2), 
Doppler boosting acts as a multiplicative factor on the light 
curve describing the flare-induced process. For a rotating jet, 
this means that Eq. 9 not only describes the oscillation of the 
average flux density of the jet, but also a modulation of the am- 
plitude AS of all short-term flares in it. For the special case of 
(1 - A)-^^" < (cr)/ (AS ), this can lead to a pattern in which short 
term variability, clearly detectable in the high states of boosting, 
becomes insignificant in the low states. For a more detailed treat- 
ment of rotating patterns with superimposed inhomogenities we 
refer to Ostorero et al. (2004), and references therein. 
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Fig. A.l. WMAP/Planck light curves for sources with at least one good timeline (i.e., four or more > 3cr data points in seven years 
of WMAP observation) in our sample. For each good timeline, all > 2cr data points are plotted, with dash-dot lines indicating the 3cr 
level. The WMAP flux densities are shown by filled symbols and the Planck flux densities are shown by open symbols of the same 
shape. 



